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EXPANDED ABSTRACT 

Two algorithms for obtaining static contact solutions are described in this 
presentation. Although they were derived for contact problems involving specific 
structures (a tire and a solid rubber cylinder), they are sufficiently general to 
be applied to other shell-of -revolution and solid-body contact problems. 

The shell-of -revolution contact algorithm is a method of obtaining a point 
load influence coefficient matrix for the portion of shell surface that is ex- 
pected to carry a contact load. If the shell is sufficiently linear with respect to 
contact loading, a single influence coefficient matrix can be used to obtain a good 
approximation of the contact pressure distribution. Otherwise, the matrix will be 
updated to reflect nonlinear load-deflection behavior. 

The solid-body contact algorithm utilizes a Lagrange multiplier to include the 
contact constraint in a potential energy functional. The solution is found by 
applying the principle of minimum potential energy. The Lagrange multiplier is 
identified as the contact load resultant for a specific deflection. 

At present, only frictionless contact solutions have been obtained with these 
algorithms. A sliding tread element has been developed to calculate friction shear 
force in the contact region of the rolling shell-of-revolution tire model. This 
element allows a relatively general, non-Coulomb, friction law to be specified for 
the contact interface. It has the added advantage of allowing friction to be cal- 
culated in the continuous interface and, when coupled with the solid-body contact 
algorithm, will permit analytic investigation of various continuum friction theories 
that have been proposed. 

The outline of future directions for the development of contact solution 
algorithms is: 


I. SHELL-OF-REVOLUTION CONTACT ALGORITHM 

II. SOLID-BODY CONTACT ALGORITHM 

III. STATIC AND ROLLING CONTACT FRICTION 

IV. FUTURE DIRECTIONS FOR CONTACT SOLUTION ALGORITHMS 
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I. SHELL-OF-REVOLUTION CONTACT ALGORITHM 


A shell whose geometry and material properties are axisymmetric can be eco- 
nomically modeled by shell-of-revolution finite elements. The SAMMSOR/SNASOR pro- 
grams (refs. 1,2), for example, permit nonlinear behavior of orthotropic shells of 
revolution to be calculated, including response to nonaxisymmetric loads. The 
algorithm described here was developed to calculate the shell deflection in response 
to a nodal point load, utilizing the calculated response to a sequence of harmoni- 
cally varying ring loads on the node. The point load solution is then used to con- 
struct an influence coefficient matrix, from which the shell contact solution is 
obtained . 



FINITE-ELEMENT TIRE MODEL 


The tire is modeled here by an assembly of axisymmetric shell elements con- 
nected to form a meridian of arbitrary curvature and following the carcass mid- 
surface. The elements are homogeneous orthotropic, with moduli determined by the 
ply structure of a particular tire. Details of this model are given in reference 3. 
If the deformation is symmetric about the wheel plane only one-half of the meridian 
is modeled, as shown in figure 1. The finite elements are /joined at nodal circles, 
referred to here as nodes. Node 12 in figure 1 is located at the tire bead and is 
given in built-in end condition. 
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Figure 1 


293 


SINGLE HARMONIC RING LOADS 


The finite-element tire model will respond to single harmonic ring loads on 
the nodal circles in addition to a uniform inflation pressure load. An approxi- 
mately linear ring load-deflection response is obtained when an individual ring load 
is applied to any node of the pressurized tire model. An example ring load- 
deflection calculation for a passenger tire model is shown in figure 2. A harmonic 
sequence of stiffness matrices is obtained by applying a sequence of single harmonic 
ring loads to each of the nodes that may be in the tire-pavement contact region. 



CROWN LOAD-DEFLECTION DATA CALCULATED WITH A UNIFORM RING LOAD APPLIED TO THE CROWN NODE 



SINGLE HARMONIC RING LOADS APPLIED TO A FINITE-ELEMENT NODE 


Figure 2 
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TRANSFER FUNCTION DEFINITION 


As a consequence of the linearity of the ring load-deflection response, the 
application of a single harmonic ring load produces a displacement field that 
varies circumferentially in the same harmonic as the applied ring load. The defini- 
tion of the transfer function T n as the ratio of the output and input amplitudes 
is given below (ref. 4). Since each node responds differently, a transfer function 
matrix T ik|n usec * to store the stiffness information generated by the ring loads. 

The partitions of this matrix are determined by the direction of the ring load. (Fig. 3) . 
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POINT LOAD VECTOR {p} AND THE DISCRETE FOURIER TRANSFORM (DFT) 

This application of the discrete Fourier transform uses an even number of 
points (N) , equally spaced around the circumference. The example shown in figure 
4 uses N=8 points. A unit load is applied at any point, say point 0. The DFT 
of the load vector yields a set of N coefficients, G . , which are approximate values 
of the coefficients of the conventional Fourier series defined on the continuous 
interval 0 £ 0 <_ 2 tt and representing the unit point load. The point load is 
applied, sequentially, in the radial, axial, and circumferential directions. 


INFLUENCE COEFFICIENT GENERAT ION 
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Figure 4 
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INVERSE DISCRETE FOURIER TRANSFORM (IDFT) AND THE INFLUENCE COEFFICIENTS 

Having the unit point load represented by a conventional Fourier series, whose 
coefficients a n are approximately given by the DFT coefficients, the transfer 
functions Tj_k| n are applied, on each harmonic, to obtain the coefficients b of the 
Fourier series representing the response of the nodal circle to the unit point 
load. The inverse discrete Fourier transform is then used to evaluate the dis~ 
placements, u , at the N points. These displacements are the elements of the 
influence coefficient matrix [A^ ik ( F± g- 5 )* 


INPUT SERIES COEFFICIENTS a = G = l 

n n N 

OUTPUT SERIES COEFFICIENTS b = a T.. , = l T.. , 
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Figure 5 
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INFLUENCE COEFFICIENT MATRIX 


The influence coefficient matrix relates the radial, axial, and circumferen- 
tial components of the displacement of points on the tire surface to the radial, 
axial, and circumferential components of load at these points. The radial response 
partition, shown in figure 6, is used to obtain a solution for frictionless con- 
tact, in which the axial and circumferential force components are known to be zero. 
The matrix here covers 3 points on each of 5 nodes. The point separation with this 
matrix is 11.25 degrees. 


P k£ = load at point l on node k 
d.j * deflection of point j on node i 


d 11 


Ann 



d 21 


A2111 A 212 i 


P 21 

d 31 


A3 111 A 3121 A 3131 { d ij} = [ A ijla] { P kt} 


P 31 

d 41 


A4III A4121 A4 131 A 4 141 


P 41 

d 51 


Asm A5J21 As 131 A 5141 A 5151 


P 51 


V 



> 

d 12 


Al 211 A 12 21 A 123 i A 124 i A 12 51 A 12 12 


P 12 

d 22 


A2211 A2221 A 2231 A 224 i A 2 251 A2212 A2222 


P 22 

d 32 


A321I A 3 22l A 323 j A 3241 A 32 51 A 32 12 A 322 2 A 323 2 


P 32 

d 42 


A42II ^4221 A 4231 A4241 A425I A 4212 A 42 22 A 4232 A4242 


P 42 

d S 2 


A52II A5221 A S231 A5241 A525J A5252 A5222 A 5 23 2 A5242 A5252 


P 52 


Figure 6 


298 



TOROIDAL SHELL CONTACT SCHEMATIC 


A cylindrical coordinate system is used to locate points on the toroidal 
surface. The coordinates r, 0, and z indicate the radial, circumferential, and 
axial directions, respectively. The tire equator lies in the r-0 plane (wheel 
plane) and a tire meridian is in an r-z plane. 

After the inflation solution has been obtained, the tire model is deflected 
against a frictionless, flat surface. The contacting surface is perpendicular 
to the wheel plane and positioned at the specified loaded radius R^, as shown 
in figure 7. The vertical load and the contact pressure distribution are unknown, 

a priori. 



Figure 7 



DEFLECTED MERIDIAN 


The deflected shape of the meridian passing through the center of contact is 
shown in figure 8. This shape is calculated by the finite-element tire model for 
the specified tire deflection of one inch. The tire load that will produce a one 
inch tire deflection is calculated to be 10,590 lb. Figure 8 also shows the 
meridian prior to inflation and the calculated shape of the meridian of the inflated 
tire, prior to contact loading. These finite-element meridians follow the carcass 
midsurface, as indicated in figure 1. Geometric and material property data on the 
Space Shuttle nose gear tire were used for the calculated results shown in figures 
8, 9, and 10. 
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Figure 8 
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CONTACT PRESSURE DISTRIBUTIONS 


The static contact pressure values (psi) calculated for two different loads 
on the Shuttle nose gear tire are shown in figure 9. The number of finite-element 
points in the contact region increases as the tire load increases. A rough esti- 
mate of the contact boundary is obtained by extrapolation of the pressure distri- 
bution. Integration of the pressure distribution gives the tire load. 


NOSE GEAR TIRE CONTACT PRESSURE DISTRIBUTIONS 
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TIRE LOAD VERSUS TIRE DEFLECTION 

An important test of a tire model is its ability to calculate a static 
load-deflection curve. Figure 10 compares the load-deflection curve calculated 
for the Shuttle nose gear tire with measured data for a similiar aircraft tire. 
Although these are both 32 x 8.8 Type VII tires, constructional details can 
alter the load-deflection curve (and many other aspects of tire behavior) . The 
cord used in the test tire is unknown and may be quite different from the nylon 
cord in the Space Shuttle tire. 


Load-Deflection Curves For 32 * 8.8 Type VII Aircraft Tires 



Figure 10 



II. SOLID- BODY CONTACT ALGORITHM 


A solid body with a cylindrical surface is often loaded in contact against a 
rigid surface. The contact load may revolve around the body, as in the case of 
a roller or a solid tire, or may remain stationary if the cylinder is used as a 
support cushion. 

Interfacial friction, present in all contact problems, is currently an active 
field of research. The study of frictional behavior is facilitated if the contact 
region is relatively large, as is produced when the body is highly deformable. 

This makes it easier to calculate distributions of normal pressure and tangential 
motion (slip) in the interface. In the case of rubber contact, the behavior 
deviates sufficiently from the Coulomb friction law that other, more physically 
realistic laws, can be easily tested. Since friction is a microscopic phenomenon, 
a contact solution giving continuous distributions of interfacial pressure and 
slip is desirable for analytic purposes. The contact algorithm described here 
provides a continuum solution for frictionless contact, the first step toward 
analysis of friction in the continuous contact interface. 


A PLANE STRAIN CONTACT PROBLEM 


An elastic semicylinder of radius R is bonded to a fixed surface, figure 
11(a). A contact load is applied by a rigid plate that deflects the semicylinder 
as shown in figure 11(b). 

The problem is formulated in terms of cylindrical material coordinates 
(r,0,z) which identify points in the undeformed body, B Q . A point P G in B Q is 
located by Cartesian coordinates x x and x 2 axes shown in figure 11(a) and 

x = r cos 0 x 0 = r sin 0 x = z 

1 2 3 

The contact load is assumed to produce a plane strain deformation. Point P Q 
moves to position P in the deformed body, B. Point P is located by the Cartesian 
coordinates y . With plane strain, y x = y x (r,0), y 2 = y 2 (r,0), and y 3 = X 3 z 
where X 3 = 1 is a specified constant extension ratio. 




Figure 11 
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GEOMETRIC DESCRIPTION 


The metric tensors g^ and G^ j , given below, completely describe the elastic 

semicylinder before and after deformation. Since y 3 is known a priori, the 

problem is solved by finding the functions y 1 (r 9 Q) and y 2 (r,0) which determine 

^± 2 * The displacement field is not utilized in this formulation but it can, of 

course, be found when x. and y. are known. (Fig. 12). 
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MATERIAL DESCRIPTION 


The material is assumed to be hyperelastic so that its constitutive properties 
are contained in a strain energy density, W. Isotropy is also assumed. For 
plane strain of an isotropic material, the strain energy is known to be a function 
of only the first and third strain invarients, Ij and I3 (ref. 5) 

W = WCl^Ig) 

For general deformation, the strain invarients are given by 


V g 1J G 


1J 


I 3 - G /g 


where g = det [g^.] = r for the semicylinder and G = det [G_]. 

When the material is also assumed to be incompressible (I 3 = 1) the constitutive 
behavior is not completely determined by the strain energy density. Hydrostatic 
pressure becomes an additional unknown, which can be determined as a Lagrange 
multiplier (ref. 6 ). This difficulty is avoided if a compressible material model 
is used. 


The material description selected for the contact problem solved here is the 
compressible neo-Hookean model developed for continuum rubber by Blatz and Ko 
(ref. 7). The Blatz-Ko model may be expressed as 


x / 3 i 

w (ii »i 3 ) = Ui - 3i 3 ) + £k| 


k-l k-l 




where y is the classical shear modulus, K is the bulk modulus, and k is a parameter 
related to atomic repulsion. When I 3 = 1, the Blatz-Ko model reduces to the 
neo-Hookean model for incompressible material. For small strains it reduces to 
the energy density giving Hooke’s law for compressible isotropic material. 
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CONTACT CONSTRAINT 


The assumption of frictionless contact with a rigid surface, positioned per- 
pendicular to the y 2 axis (see fig. 11(b)), implies a geometric constraint only on 
the solution function y^(r,0). Within the contact region, whose extent is not 
known a priori the deformed surface is flat and it is known that y 2 = 
where R £ is the specified location of the contact surface. Outside the contact 
region, and in the interior of B, the solution must satisfy y 2 < R£. This in- 
equality constraint on y is converted to an equality constraint by introducing a 
new function s(r,0), defined by the following equation 

2 

y^ + s = R^ (constraint equation) 

which is valid everywhere on the boundary and in the interior of B. The function 
s(r,0), called a slack variable, has been used previously in optimization problems 
with an inequality constraint. Reference 8 gives several examples of the use of 
slack variables. 


The contact problem is solved by minimizing the strain energy in B, subject 
to the constraint equation given above. Since plane strain is assumed, the energy 
is uniform along the axis of the semicylinder. Using symmetry, integration of the 
energy density is taken over one-half of the r-0 plane contained in B^. 


The constraint equation is brought into the energy density functional by 
means of a Lagrange multiplier function A(r,0). The contact problem is then 
governed by the following functional 


where 


Uyj. y 2 , s, X) 



e, y. 


y 2 , s, X) drd0 


F = rV ^ y a,^ + X + s2 


Although rW is positive definite, F is not positive definite due to the addi- 
tion of the constraint. Therefore, I may only be regarded as being made sta- 
tionary instead of minimized by equilibrium solution functions. The integral 
of W, however, is minimized by the equilibrium solution and this is used as a 
check during the solution finding process. 


Through additional analysis (J. T. Tielking, Texas A and M University, 
unpublished data) the Lagrange multiplier function is shown to be an unknown 
constant, identified as the resultant load in the contact region. The constraint 
condition may then be removed from the integral, the slack variable is no longer 
needed, and the contact problem is now governed by 

Ky a = JJ rW(y a drd6 + X (y 2 (R,tt/2)-R^) 

The solution is obtained by finding y a (r,0) and the constant X which make 
I(y a » X) stationary. (Fig. 13). 
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SOLUTION FUNCTIONS 


The numerical solution is obtained by application of the principle of 
stationary potential energy (ref. 6), using the functional I(ya, X) in which X is 

an unknown constant. The solution functions y are taken as two-dimensional 

0 1 

finite series 

M N 

yj (r,0) = Xl (r,0) + ^2 2 a ij y i ±j (r,0) 

i-1 j=l 

M N 

y z (r,e> = x 2 (r,6) + ^2 b ij y 2 lj (r,0) 

i=l 3=1 

where a., and b.. are unknown coefficients. The symmetry and geometric boundary 
conditions evident in figure 11 are met by 

y ^ (r,0) = r 1 sin (2j0) 

ij i 

y (r,0) = r sin [(2j- 1)0] 

2 

The above functions allow the energy density to be integrated, thereby re- 
ducing the functional I to an algebraic function of the 2 x M x N +1 unknown 

constants a.., b.., and A 
ij iJ 

l - l («y ; by ; X) 

The functional I is made stationary by the constants obtained from the following 
set of simultaneous nonlinear equations: 

9 1 t 

"KZ = 0 and — — = 0 for i = 1,2, . . . ,M and j = 1,2, . . ., N 
U 3by 

|i = y 2 (R, 71/2 ) - = 0 

This system is solved in an iterative manner by the Newton-Raphson method. Using 
the starting values a.. = b., = X = 0, five or six iterations (which give successive 
corrections to these constant's) are usually sufficient. The iterations are con- 
tinued until the corrections appear to have negligible effect on the solution 
functions y a (r,8). The energy density is evaluated after each iteration to 
check for minimization. 
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DEFORMATION SOLUTION 


Numerical results have been obtained using material constants y = 75 psi 
(shear modulus), K = 475,000 psi (bulk modulus), and k = 13,3 in the Blatz-Ko 
model. The values of K and k are taken from reference 7 where they are shown to 
give a good fit to hydrostatic compression data on Butyl tread rubber (polyiso- 
butylene). The shear modulus is believed to be a realistic estimate, based on 
Treloar's statement (ref. 9) that the shear modulus of rubber is lower than the 
bulk modulus by a factor of about 10^. 

The computer-generated drawing below shows coordinate circles and radii 
before deformation (x x and x 2 , dashed lines) and the deformed configuration 
(y^ and y 2 , solid lines) of these circles and radii. The deformation is pro- 
duced by a 10-percent deflection of the contacting surface (shown dashed). The 
deformation solution, y x (r,0) and y 2 (r,0), is obtained in a 16-term series 
for each function; the strain energy is minimized by the coefficients a^. and 
b^j for i=l,2,3,4 and j=l,2,3,4, found after six Newton-Raphson iterations. 

This computation took 30 seconds of CPU time on a mainframe computer (Amdahl 
470/V8). The Lagrange multiplier obtained in this solution is A=93.1, inter- 
preted as a 93.1 lb load needed for a 10-percent deflection if the semi cylinder 
extends 1 inch in the z-direction. (Fig. 14) . 
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III. STATIC AND ROLLING CONTACT FRICTION 


Static Contact Friction. A body is brought into static contact by motion 
perpendicular to the contact plane. During this motion, the contact boundary 
expands until the resultant of the normal contact pressure reaches equilibrium 
with the external load applied to the body. As the contact region is formed, shear 
forces are generated by tangential motion of contacting surface points. These 
shear forces are frictional and transient, reaching equilibrium levels when the 
body itself comes into equilibrium. Although the body may be assumed elastic, and 
thus conservative, the frictional shear forces are not conservative. 

The formidable problem of calculating a static contact solution including the 
effect of friction is alleviated somewhat by assuming Coulomb’s law of friction 
is valid in the contact region. An algorithm for including Coulomb friction in 
a static contact problem has been developed by Rothert et al. (ref. 10). Although 
Coulomb friction may be taken for an approximate analysis of frictional contact, 
mathematical and physical uncertainties arise when it is assumed. Nonclassical 
friction laws have been proposed by Oden and his coworkers (e.g., ref. 11). Al- 
though these appear to have been developed mainly for metals, they may be applic- 
able to more deformable material such as rubber. 

Rotting Contact Friction . This discussion is limited to steady rolling 
under a constant load. Neglecting hysteretic effects in the body, the power in- 
put to maintain steady rolling is balanced by the work rate of the friction forces 
in the contact region. In steady rolling, a contact region, whose boundary is 
fixed by the load, is continuously generated. The normal pressure and sliding 
velocity at a given location within the contact boundary do not change with time 
so steady-state frictional behavior is maintained. The rolling contact problem 
with friction is therefore much easier to analyze and provides a mechanism for 
the study of nonclassical friction theories. An algorithm for calculating friction 
in the contact region of a tire rolling at constant velocity will be described 
next. 
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SLIDING TREAD MODEL 


This is a tread element model developed to convert a frictionless rolling 
contact solution into a solution for rolling contact with friction. The element 
passes through the contact region with the velocity found for frictionless contact. 
This is termed the carcass velocity, V_ , whose distribution is symmetric about the 
center of the footprint as sketched below 




For a free-rolling tire, the amplitude of the footprint sliding velocity is 
very small. At 60 mph (1056 ips) , the peak V c is calculated (by the author) to 
be about 50 ips in a frictionless footprint. 

In free-rolling, the normal contact pressure distribution, p, is essentially 
unchanged by friction. The sliding velocity distribution, however, is signifi- 
cally altered in an interactive manner. A hysteretic theory of tire-pavement 
friction proposed by Schapery (ref. 12) gives the dependence of the friction 
force, F s , on the actual sliding velocity, V s , and normal pressure, p, at a 
point in the contact region. This is expressed as 

I 3 l b 

F = - B sgn(V ) j V I x P (hysteretic theory) 

s s s 

where a, b, and B are material friction properties. The sliding tread element 
model, shown in figure 15(a), is viscoelastic with stiffness and damping parameters 
K and c. Sliding friction, F , causes the element to deform, thereby influencing 
the sliding velocity V g . The following nonlinear differential equation is 
integrated to calculate V g . 

dV K(V -V ) + cV 

s __ c s c__ 

dt ' dF s 



In this equation, dF s /dV g is the rate of change of sliding friction with sliding 
velocity. This can be obtained by differentiating the hysteretic theory given 
above or measured experimentally. Footprint transit time, t, is taken as the 
independent variable. The time is equivalent to location in the footprint for 
steady rolling. 
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IV. FUTURE DIRECTIONS FOR CONTACT SOLUTION ALGORITHMS 


As in other areas of solid mechanics, future research on contact problems 
will be directed towards obtaining solutions valid for large deformations. The 
finite-element method seems particularly well suited for application to contact 
problems and special elements have already been developed for this purpose. Con- 
tinuum mechanics research on contact problems should not be neglected, however. 

A large-deformation contact solution in terms of continuous functions will prove 
valuable in the analysis of contact with friction and the assessment of friction 
laws now being proposed for deformable bodies. 

A true contact problem is one in which the contact boundary and interfacial 
pressure distributions are unknown a priori. At present, it appears that such 
problems will be displacement prescribed: Deflection of the body toward the 

contact surface is specified and integration of the calculated normal component 
of the contact pressure gives the resultant load. Some effort should be directed 
towards a load-specified contact problem, perhaps utilizing the principle of 
stationary complementary energy to calculate the interfacial pressure distribu- 
tion subject to the prescribed load constraint. Validated solutions for fric- 
tionless contact are essential prior to including the effect of friction on the 
contact solution. 

In the analysis of friction, it seems that the study of rolling contact as 
a steady-state problem has much to offer. As friction is an interactive phe- 
nomenon, at least in regard to sliding velocity, sophisticated algorithms are 
needed to generate the frictional contact solution from the solution for fric- 
tionless contact (which will undoubtably be the starting point). 

The following schematic, figure 16, outlines a progression of research on 
contact problems. Linear sliding contact is excluded from the outline as this 
is usually a transient situation leading to accelerated wear and abrasion. 

Much more will be gained by research focused on rolling contact which will, in 
any case, include sliding. 
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PROGRESSION OF RESEARCH ON FRICTIONAL CONTACT PROBLEMS 



Figure 16 
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